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ABSTRACT 

The  following  report  defines  the  scope,  provides  a  narration  of 
the  technical  accomplishments  and  current  status  after  one  year  of 
operation,  and  lists  the  work  remaining  to  be  done  under  U.S.  Bureau 
of  Mines  Contract  No.  H0210022  with  the  University  of  California, 
Berkeley,  on  a  project  entitled,  "Wave  Propagation  in  Anisotropic 
Rocks,"  The  project  has  been  divided  into  throe  major  phases,  namely 
(a)  an  experimental  examination  of  the  wave  processes  in  Yule  marble 
blocks,  (b)  a  theoretical  study  of  the  pulse  propagation  produced  by 
a  concentrated  transient  source  normal  to  the  surface  of  a  homogeneous, 
elastic  transversely  isotropic  semi-infinite  solid  using  (1)  an  integral 
transform  technique,  and  (2)  a  finite  element  method,  and  (c)  a 
determination  of  the  geometric  and  static  and  dynamic  mechanical 
properties  of  tne  Yule  marble,  including  fracture  characteristics. 

In  the  experimental  wave  propagation  phase,  the  major  accomplish¬ 
ments  achieved  to  date  include  the  development  of  a  crystal  transducer 
package  with  a  laterally  unconstrained  crystal  employed  as  a  sensing 
element,  the  calibration  of  both  the  crystal  and  the  entire  package, 
the  development  of  installation  techniques  for  the  embedment  of  the 
transducer  inside  cores  drilleu  in  rock  bars,  and  some  progress  in  the 
development  of  a  suitable  grouting  material  to  fill  the  core  holes  after 
installation  of  the  transducer  unit  with  a  minimum  of  dynamic  mismatch. 
The  integral  transform  technique  has  been  developed  and  partially  tested 
to  obtain  the  displacement  field  in  the  transversely  isotropic  solid 
with  an  axis  located  in  the  free  surface  under  a  Heaviside  input,  This 
process  constructs  a  three-dimensional  solution  from  the  integration  of 
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a  series  of  two-dimensional  problems  associated  with  line  loads  on  the 
surface  of  the  half-space.  It  employs  a  Cagniard-de  Hoop  transforma¬ 
tion  which  simplifies  the  inversion  process  and  leads  to  a  physical 
interpretation  of  the  transient  wave  process  in  terms  of  well-established 
concepts  of  wave  and  slowness  surfaces  that  have  been  employed  in  the 
field  of  crystal  acoustics. 

A  finite  element  program  has  been  written  for  the  three-dimensional 
problem  cited  above  and  is  currently  in  the  debugging  and  test  phase. 

It  has  been  specialized  to  the  simpler  case  of  inotropic  behavior  where 
other  solutions  for  checking  purposes  exist.  Results  obtained  tc  date  include 
the  uniaxial  wave  process  in  a  rod  and  the  surface  motion  of  a  half¬ 
space.  They  are  very  encouraging  in  terms  of  providing  good  correspon¬ 
dence  with  these  other  solutions,  and  indicate  the  immediate  application 
of  the  program  to  the  subject  problem  after  some  additional  investiga¬ 
tion  of  the  numerical  integration  schemes. 

Crystallographic  techniques  have  been  developed  for  the  location 

of  the  axis  of  elastic  symmetry  of  the  Yule  marble  specimens.  Static 

compressive  tests,  some  with  repeated  loading,  have  been  conducted  on 

samples  of  the  material,  indicating  significant  non-linearities  and 

the  presence  of  hysteresis,  A  technique  has  been  developed  for  the 

generation  of  constant  strain  rates  in  the  intermediate  range  of  10 
-X 

to  100  sec  utilizing  an  adaptation  of  a  Hopkinson-bar  procedure. 

This  work  will  be  continued  and,  during  the  coming  year,  the 
various  phases  of  the  program  will  be  intimately  correlated.  The  pro¬ 
gram  is  on  schedule  and  no  major  difficulties  are  anticipated  in  the 
completion  of  the  work  according  to  the  original  schedule. 
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PREPACE 

This  is  the  first  annual  report  concerned  with  the  technical 
activities  performed  under  Contract  No.  R021QG22  of  the  U.S.  Bureau  of 
Mines  with  the  University  of  California  on  the  subject  entitled  "Wave 
ropagation  in  Anisotropic  Rocks."  It  covers  the  period  from  February 
5,  1971  to  February  4,  1972  and  will  be  followed  by  additional  reports 
involving  extensions  of  the  present  contract.  Ihe  scope  of  the  program 
is  detailed  in  the  proposal  identified  as  UCB-Eng  3154,  dated  19  August, 
1970  and  submitted  on  behalf  of  the  University  of  California,  Berkeley, 
by  W.  Goldsmith  as  Principal  Investigator,  and  is  also  spelled  out  in 
Article  I  of  the  subject  contract.  It  was  envisioned  initially  that 
the  completion  of  the  entire  investigation  outlined  in  these  documents 
would  require  a  period  somewhat  in  excess  of  two  years;  consequently, 
the  present  document  comprises  an  interim  rather  than  a  terminal  report 
of  the  program.  The  time  schedule  originally  indicated  for  the  various 
phases  of  the  work  has  been  closely  observed. 

The  scope  of  the  activities  was  subdivided  into  several  categories 
as  follows;  (a)  An  experimental  examination  of  the  wave  processes 
produced  by  impacts  of  spheres  on  the  surface  of  large  blocks  of  Yule 
marble,  involving  the  measurement  of  surface  and  body  waves  by  means 
of  strain  gages,  accelerometers,  and  transducer  packages  embedded  in  the 
solid;  (b)  A  theoretical  analysis  of  the  pulses  produced  in  the  block, 
modelled  as  a  transversely  anisotropic  half-space,  using  integral 
transform  and  other  mathematical  techniques,  and  comparison  of  the  pre¬ 
dicted  response  of  this  model  with  experimental  data;  and  (c)  The  exam¬ 
ination  of  the  mechanical  properties  of  the  material  Yule  marble. 
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including  th«  effect  of  rate  of  loading  on  fracture  properties.  Each 
of  four  graduate  students  was  assigned  to  a  separate  phase  of  the  program: 
Mr.  K.  Krlshnaaoorthy  devoted  his  efforts  to  part  (a);  Mr.  S.  h.  Suh 
was  concerned  with  the  development  and  numerical  evaluation  of  the 
analysis  employing  integral  transform  techniques;  Mr.  M.  Katona 
occupied  himself  with  the  development,  validation  and  application  o'  a 
finite-element  solution  oi  the  problem;  and  Mr.  S.  Howe  was  associated 
with  the  phase  involving  the  determination  of  the  mechanical  properties 
of  the  rock  including  fracture.  These  students  are  utilizing  the  re¬ 
sults  of  their  efforts  on  this  project  in  whole  or  in  part  for  their 
doctoral  dissertations  or  masters  theses.  Professor  W.  Goldsmith, 

Division  of  Applied  Mechanics,  and  Professor  J.  L.  Sackman,  Division  of 
Structural  Engineering  and  Structural  Mechanics,  jointly  supervised  the 
entire  program  since  its  .nception  and  are  responsible  for  its  conduct. 
Technical  assistance  in  certain  portions  of  the  Fortran  programming  was 
rendered  by  Mr.  R.  Kenner,  and  the  technical  staff  of  the  Departments 
of  Mechanical  Engineering,  Civil  Engineering  and  Geology  and  Geophysics 
assist  d  in  the  operation  of  the  experimental  investigation. 

There  have  been  no  major  technical  difficulties  in  the  pursuit  of 
the  original  goals,  and  the  progress  made  is  encouraging  to  the  point  of 
predicting  that  these  will  be  achieved  within  the  original  time  span 
estimated.  Some  changes  in  emphasis  or  techniques  may  be  required, 
particularly  in  the  methodology  of  determining  the  strain-rate  effects 
on  the  fracture  process.  Furthermore,  an  additional  effort  may  be  made 
to  provide  a  sound  mechanical  model  for  the  Yule  marble.  The  program 
has  already  Indicated  the  desirability  of  additional  investigations  in 
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the  dOMin  of  anisotropic  wave  propagation  that  will  be  detailed  in 
the  appropriate  section  of  the  report.  The  prog raw  has  exhausted  the 
funds  initially  requested  without  any  overrun. 
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I.  IKTROPPCTIOM 
a)  Htatorlcal  Background 

A  detailed  examination  of  the  analyais  of  transients  in  homogeneous, 

anisot Topically  elastic  media  has  been  undertaken  only  during  the  last 
(1-71* 

fifteen  years  ,  ind  specific  nuaerical  results  have  been  obtained 

essentially  only  for  two-diaensional  classes  of  probleats.  One  of  these 
is  the  solution  to  the  application  of  a  tine-dependent  point  load  on 
the  free  suriace  of  a  quasi-isotropic  half~space,  i.e.  a  mnterial  in 
which  there  exists  an  axis  of  elastic  symmetry  perpendicular  to  the 
tree  surface,  leading  to  an  axisymmetric  problem.  Another  result 

(7) 

has  been  recently  obtained  by  Burridge  for  a  line  load  applied  to 
the  surface  of  a  semi-inf inite  solid,  with  cubic  symmetry,  i.e.  one 
wnere  the  axis  of  elastic  symmetry  lies  in  the  free  surface;  this  is 
also  a  two-dimensional  problem,  albeit  of  a  different  nature.  In  spite 
of  the  limited  nature  of  the  presently  available  theoretical  develop¬ 
ment,  there  have  been  no  corresponding  experimental  investigations 
whose  results  could  be  compared  with  the  predictions  of  these  theories. 

As  an  initial  attempt  to  broaden  the  scope  of  knowledge  in  this 
domain,  a  closely  coordinated  tneoretical  and  experimental  program  was 
undertaken  at  the  University  of  California  by  Dr.  T,  E.  Ricketts  that 
was  concerned  with  the  surface  wave  pulse  propagation  produced  by  normal 

impact  of  a  sphere  on  a  homogeneous,  transversely  anisotropic 
(8) (9) 

half-space.  For  the  experimental  phase,  this  solid  was  modelled 

by  a  large  block  of  Yule  marble  whose  axis  of  elastic  symmetry  was 


*Numbers  in  superiors  refer  to  the  Bibliography  at  the  end  of  the  report. 
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located  in  the  free  surface  plane.  Analytical  expressions  were  ob¬ 
tained  using  a  combined  Laplace- Fourier  transform  method  which  provided 
the  formal  solution  of  the  three-dimensional  problem  in  transform 
space.  The  relations  were  derived  for  general  anisotropy  and  then 
specialized  for  the  transversely  anisotropic  case;  with  the  elastic 
symmetry  axis  in  the  free  surface  however,  a  complete  inversion  of  the 
displacement  field  into  the  real  time-space  domain  was  not  obtained. 

The  experiment  utilized  3/4"  diameter  steel  spheres  which  struck 
the  test  block  at  impact  velocities  of  about  20  ft/sec.  Semi-conductor 
strain  gages  and  specially-fabricated  quartz  crystal  accelerometers 
were  employed  to  detect  the  horizontal  and  vertical  components  of  the 
disturbance,  respectively.  The  former  were  mounted  both  along  and 
normal  to  radial  lines  drawn  from  the  impact  point  in  order  to  detect 
any  quasi-Rayleigh  components. 

The  experimental  procedure  provided  surface  wave  group  velocity, 
group  slowness,  amplitude  and  attenuation  curves  for  the  Yule  marble. 
Hie  experimental  slowness  curve  was  compared  to  the  corresponding 
analytical  phase  and  group  slowness  curves  computed  using  both  static 
and  dynamic  properties  of  this  substance  which  were  determined  i  i  this 
study.  A  comparison  of  analytical  and  experimental  results  provided 
acceptable  correlation. 

The  investigation  described  above  represented  the  first  combined 
analytical  and  experimental  program  involving  transient  response  of  ,ny 
anisotropic  elastic  solid  and  also  the  first  to  consider  a  fulli  three- 
dimensional  problem.  TTiis  pioneering  effo.-t  w  s  successful,  but  ■»! 


necessity  of  limned  sc»pe,  providing  basic  features  of  the  solution 
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applicable  to  surface  wave  phenomena.  Details  of  the  body  wave  motion 
were  not  obtained,  but  the  work  represents  a  foundation  for  future  in¬ 
vestigations  and  served  as  the  motivation  for  the  scope  of  the  current 
contract. 

b)  Scope  of  the  Research  Proposed  Under  Contract  Ho.  H0210Q22 
The  purpose  of  the  investigation  was  the  determination  of  the 
response  of  a  re.’l  rock  material  to  a  surface  impact.  To  this  end,  the 
variation  of  stresses  in  the  interior  and  on  the  surface  of  a  block  of 
Yule  marble  to  transient  loading  by  a  spherical  projectile  were  to  be 
measured  and  c  imp a red  to  a  corresponding  mathematical  model  of  the 
system. 

Two  large  blocks  of  Yule  marble  were  supplied  by  the  U. S.  Bureau 
of  Mines  for  testing  purposes.  This  material  was  selected  since  it 
represents  a  reasonably  homogeneous  continuum  or»  the  scale  of  the 
wavelengths  employed,  yet  exhibits  a  greater  degree  of  transverse 
anisotropy  than  any  other  natural  substance  available  in  the  required 
size  and  quality.  In  one  of  the  blocks,  the  axis  of  elastic  symmetry 
is  in  one  of  the  nee  surfaces,  but  not  parallel  to  the  edges; 
consequently,  the  block  can  be  used  either  as  a  transversely  isotropic 
solid  loaded  in  a  direction  parallel  to  one  of  the  non-symmetric  axes 
or  at  an  angle  to  these.  To  be  employed  in  an  experiment  involving  a 
quasl-isotropic  situation,  it  would  be  necessary  to  machine  the  block 
so  as  to  produce  a  loading  plane  perpendicular  to  the  axis  of  elastic 
symmetry. 

The  methodology  for  the  measurement  of  surface  phenomena  has  been 
satisfactorily  worked  out  by  Ricketts  utilizing  semi-conductor  strain 
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(8) 

gages  and  accelerometers'  .  The  present  objective  rt quires  the 
development  and  Installation  of  transducers  capable  of  recording  the 
response  in  various  directions  in  the  interior  of  the  block.  To  this 
effect,  it  was  necessary  to  institute  a  developmental  program  for  the 
manufacture,  calibration  and  installation  of  minute  sensor  packages  to 
be  placed  at  the  terminus  of  various  holes  of  very  small  diameter 
relative  to  the  block  size.  In  addition  to  the  actual  problems  of 
excavation,  there  is  the  crucial  problem  of  selection  of  a  suitable 
filler  material  that  anchors  the  transducer  package,  completely  fills 
the  drill  hole  and  exhibits  appropriate  acoustic  impedance  and  other 
mechanical  characteristics  compatible  with  the  parent  material.  The 
proper  functioning  of  the  externally  calibrated  transducer  unit  in  situ 
will  require  a  special  investigative  subprogram  and  constant  attention. 

In  view  of  the  limited  space  available  for  the  placement  of  the 
transducers,  their  locations  within  the  block  must  be  selected  with 
great  care  so  as  to  permit  the  recording  of  data  for  a  wide  variety  of 
external  loading  situations  and  for  different  response  quantities, 
ta’cing  full  advantage  of  the  symmetry  of  the  block.  Loading  will  occur 
by  means  of  sphere  impacts  using  either  a  drop  or  pneumatic  techniques; 
the  initial  and  terminal  velocities  will  t>e  monitored  so  as  to  permit 
control  of  the  characteristics  of  the  incident  impulse. 

The  mathematical  model  employed  in  this  study  consists  of  a  trans¬ 
versely  isotropic  homog*neous  elastic  half-space  subjected  to  a  transient 
point  load  applied  normal  to  the  free  surface  that  contains  the  axis  of 
elastic  symmetry.  The  Integral  transform  techniques  initiated  in 
References  (8)  and  (9)  were  to  be  further  developed  so  as  to  permit 
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numerical  evaluation  using  high-speed  digital  computing  equipment  of 
required  field  variables  at  selected  stations.  During  this  process, 
it  was  found  to  be  necessary  to  employ  a  different  transformation  than 
sketched  earlier  in  order  to  obtain  a  more  practical  solution  algorithm. 
Eventually,  the  displacement  and  stress  fields  are  to  be  generated  from 
this  analysis  for  arbitrary  loading  histories. 

It  was  deemed  necessary  to  initiate  an  independent  method  of 
analysis  involving  the  discretization  of  the  block  by  means  of  the 
finite  element  method  coupled  with  a  variational  principle  from  which 
the  equations  of  motion  are  derived.  This  technique  will  be  used  to 
validate  the  results  obtained  from  the  Integral  transform  procedure  and 
hafl  a  potential  for  broader  application.*  in  the  theoretical  investiga¬ 
tion  of  more  general  wave  ,  repagation  processes  such  as  might  be 
encountered  In  the  field.  These  complications  might  include  the  effects 
of  total  anisotropy,  inhomogeneity,  the  presence  of  boundaries,  more 
general  loading  conditions,  and  possibly  inelastic  behavior.  Ibis 
effort  has  paid  unanticipated  dividends  by  both  providing  a  degree  of 
confidence  in  the  transform  approach  upon  successful  correlation  of 
corresponding  numerical  results,  and  by  indicating  the  direction  to  be 
followed  analytically  in  the  pursuit  of  the  solution  of  more  complicated 
phenomena  beyond  the  s-one  of  the  present  contract. 

In  order  to  predict  with  confidence  the  diverse  phenomena  expected 
to  occur  upon  loading  of  this  material,  it  is  crucial  that  a  full 
knowledge  of  the  complete  mechanical  behavior  and  its  dependence  upon 
loading  rate  as  well  as  the  orientation  of  the  axes  of  elastic  symmetry 
within  the  two  actual  blocks  be  precisely  delineated.  Some  information 
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on  the  values  of  the  static  and  dynamic  elastic  constants  as  well  as 
the  orientation  of  the  principal  axes  for  one  of  the  blocks  was  indi¬ 
cated  in  Reference  (8);  however,  the  precision  of  the  numerical  procedures 
to  be  employed  in  the  analysis  requires  a  verification  and  refinement 
of  these  values  as  well  as  a  more  extensive  investigation  of  the  stress, 
strain,  strain-rate  behavior  of  the  material,  composing  both  blocks.  In 
addition,  the  effect  of  loading  rate  on  the  fracture  characteristics 
of  the  material  was  to  be  determined  so  as  to  determine  the  limits  of 
applicability  of  the  analysis  detailed  above.  The  data  will  be  employed 
to  construct  a  physical  model  of  the  behavior  pattern  of  the  Yule 
marble  that  might  be  useful  in  assessing  the  deficiencies  of  the  homo¬ 
geneous  linearly  elastic  comportment  assumed  in  the  mathematical 
treatment  of  the  impact  phenomenon,  as  well  as  being  of  considerable 
value  in  its  own  right. 

Stress-strain  curves  for  the  marble  have  been  obtained  at  various 
low  rates  of  loading  and  for  repeated  loading  below  the  failure  point. 
Methods  of  achieving  higher  loading  rates  have  been  investigated  and 
will  be  applied,  A  technique  for  the  determination  of  fracture 
characteristics  of  bone^0^  will  be  utilized  for  the  corresponding  in¬ 
vestigation  of  the  marble,  supplemented  by  suitable  visual  techniques. 

A  significant  portion  of  the  total  effort  was  devoted  to  the  evolu¬ 
tion  and  development  of  new  techniques  necessary  to  achieve  both  the 
experimental  and  numerical  objectives  of  the  program.  This  has  involved 
the  selection,  design,  construction  and  calibration  of  the  transducer 
package,  the  methodology  of  drilling  and  the  fabrication  and  injection 
of  filler  material  for  the  holes.  Further,  this  includes  the  writing 
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and  debugging  of  the  computer  programs  and  their  specialization  to 
problems  with  known  solutions  for  the  purpose  of  .aiidation.  Investiga¬ 
tion  of  possible  methods  of  determining  the  principal  axes  of  the 
blocks  with  sufficient  accuracy,  yet  reasonable  in  terms  of  cost  and 
time  expenditure,  1.  ive  received  special  attention.  These  and  other 
required  preliminary  efforts  do  not  manifest  themselves  in  immediate 
concrete  results,  but  are  vitally  necessary  to  the  successful  execution 
of  the  project. 

c)  Current  Status  of  the  Project 

After  one  year  of  operation,  the  accomplishments  of  the  project 
correspond  quite  closely  to  those  antictnated  in  the  program  schedule 
of  the  original  proposal.  The  transducer  development,  budgeted  at  six 
months,  has  been  completed  and  a  calibration  procedure  standardized. 
Installation  techniques,  Including  coring  methods  and  investigation  of 
filler  materials,  are  still  in  the  process  of  examination,  alchough  the 
majority  of  the  effort  for  this  phase  has  been  expended.  When  this 
task  is  completed,  expected  to  occur  in  another  two  to  three  months, 
locations  of  the  transducers  in  the  two  Yule  marble  blocks  together 
with  access  holes  will  be  selected  and  in  situ  calibration  and  testing 
will  be  initiated. 

The  solution  of  the  boundary  value  problem  by  integral  transform 
methods  has  been  reformulated  in  terms  of  variables  different  from  those 
indicated  in  Reference  (9).  This  has  permitted  the  delineation  of  the 
displacement  field  in  terms  of  quantities  derivable  from  the  roots  of  a 
sextic  algebraic  equation  whose  coefficients  depend  upon  the  elastic 


constants  and  on  the  spatial  transform  parameters.  A  program  has  been 
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written  to  carry  out  the  necessary  numerical  computations  and  test  ■ 

cases  are  currently  under  scutiny.  Consequently,  this  aspect  of  the  . 

investigation  is  approximately  three  months  ahead  of  schedule.  An 
additional  effort  will  be  required  to  determine  an  analytical  expres¬ 
sion  for  the  strain  field  and  to  write  a  program  to  evaluate  tins 
field,  so  that  the  results  so  obtained  can  be  compared  directly  with 
strain  measurements.  This  procedure  will  probably  involve  a  large 
amount  of  computational  effort,  but  this  appears  to  be  unavoidable  since 
the  alternative  of  numerical  differentiation  is  not  a  feasible 
procedure. 

A  three-dimensional  finite  element  code  for  a  general  anisotropic 
elastic  medium  has  been  programmed  and  found  to  work  effectively  in 
predicting  the  effects  of  uniaxial  wave  propagation  in  isotropic 
materials.  The  results  for  the  three-dimensional  case  of  a  point  load 
acting  on  an  isotropic  half-space  have  been  found  to  disagree  signi¬ 
ficantly  with  those  obtained  from  the  corresponding  closed-form 
solution.  The  indications  art-  that  the  difficulties  encountered  can 
be  attributed  to  the  numerical  integration  scheme  employed  over  the 
basic  finite  element.  This  scheme  will  be  revised  and,  further,  the 
finite  element  will  be  improved.  It  is  expected  that  upon  the  resolu¬ 
tion  of  this  difficulty,  no  significant  problems  will  occur  when  the 
anisotropic  model  used  in  the  present  investigation  is  introduced  into 
the  progrn.ii. 

The  investigations  involving  the  property  determinations  of  the 
blocks  are  slightly  behind  schedule  since  this  phase  was  initiated 


almost  one-half  year  subsequent  to  the  beginning  of  the  contract. 
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However,  since  ouny  of  the  techniques  to  be  employed  here  have  pre¬ 
viously  been  utilized  extensively  in  our  laboratories,  the  overall 
development  time  for  this  portion  will  be  considerably  less  than  for 


i 


the  other  experimental  phases,  permitting  the  expenditure  of  more 
time  for  the  gathering  of  data. 

The  following  section  will  present  in  detail  the  technical  accom¬ 
plishments  during  the  first  year  of  operation  of  the  contract.  A 
subsequent  section  will  discuss  the  anticipated  efforts  during  the 


second  year  of  the  investigation. 
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II,  TECHNICAL  ACCOMPLISHMENTS 

The  achievements  attained  during  the  first  year  of  operation  of 
the  contract  are  described  in  detail  in  the  following  section.  The 
subject  aatter  is  subdivided  into  the  following  categories:  (a) 
Gxperiaental  Exaaination  of  the  Wave  Processes  in  an  Anisotropic  Rock 
Material;  (b)  Theoretical  Exaaination  of  the  Pulse  Propagation  Produced 
in  an  Anisotropic  Solid  by  a  Concentrated  Transient  Source:  (1)  Integral 
Transform  Method,  and  (2)  Finite  Element  Method;  and  (c)  Determination 
of  the  Mechanical  and  Geometrical  Properties  of  the  Yule  Marble  Blocks 
to  be  Employed  as  Half-spaces, 


a)  Experimental  Examination  of  the  Wave  Processes  in  an  Anisotropic 
Rock  Material  ~~  ‘  ““ 

In  the  present  project,  the  examination  of  body  and  nurface  waves 

will  be  conducted  using  the  two  blocks  of  Yule  marble  cited  previously, 

with  dimensions  of  l£'  X  lj '  X  2' ,  one  of  which  was  used  in  the 

(8) (9) 

investigation  by  Ricketts  ,  The  procedures  to  be  employed  also 

closely  parallel  this  earlier  work,  but  Involve  significant  and 
Innovative  extensions. 

The  specimens  will  be  loaded  on  the  free  surface  by  means  of  sphere 
impacts  either  directly  or  through  an  intermediate  loading  bar;  the 
loading  pulse  will  be  monitored  by  crystals  or  strain  gages  at  the 
contact  point  or  on  this  intermediate  bar.  The  impact  will  be  controlled 
in  order  to  limit  the  length  of  the  pulse  to  a  dimension  small  compared 
to  the  size  of  the  block  so  that  observations  of  the  entire  initial 
transmitted  pulse  at  a  number  of  stations  can  be  accomplished  without 
disturbance  due  to  waves  reflected  from  the  free  surfaces  of  the  specimen. 
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Sueh  signals  fro#  both  surface  and  Interior  detectors  will  be  photo** 
graphically  recorded  using  oscilloscopes.  Pulse  shape  variations  can 
be  accomplished  by  the  use  of  spheres  of  different  sizes,  other  shapes 
and  sizes  of  strikers,  eeployaent  of  projectiles  composed  of  different 
materials,  or  by  changing  the  shape  and/or  composition  of  the  trans¬ 
mission  devices. 

The  initial  effort  to  produce  body-wave  sensors  utilized  1/8"  X 

1/32"  circular  X-cut  quartz  crystals  that  measure  uniaxial  normal  stress. 

Ihe  first  experiments  were  performed  in  rock  bars  of  about  3/4  inch  in 

diameter  subjected  to  central  longitudina*  impact.  This  involved  the 

attachment  of  lead  wires  to  the  crystals  using  Electrobond  #2016 

adhesive  manufactured  by  Adhesive  Engineering  Company.  A  bar  about  2 

feet  long  was  cut  into  two  by  means  of  a  special  rock  saw  and  a  hole 

3/16"  in  diameter  and  3/16"  deep  was  drilled  on  the  centerline  of  one 

of  the  segments.  A  diametral  groove  was  sliced  across  the  hole  to 

accommodate  the  lead  wires,  as  shown  in  Figure  1.  The  hole  wt.s  now 

filled  with  a  mixture  of  EPY-150  epoxy  cement  and  rock  powder,  and  the 

crystal  was  embedded  in  this  composite.  The  other  segment  was  then 

bonded  to  the  first  to  produce  a  single  test  specimen  with  an  internal 

(8) 

transducer,  employing  a  special  alignment  jig  .  A  pair  of  SP.-4, 
FAE-12-12  S9L  foil  gages  were  mounted  on  the  opposite  ends  of  a  bar 
diameter  at  a  position  about  3  inches  from  the  crystal  position  by  means 
of  EPY-150  cement  manufactured  by  Baldwin-Lima-Hami11  ton  Electronics, 

Inc.,  as  shown  in  Figure  2.  These  gages  were  incorporated  in  a 
potent iometric  circuit  and  w<jre  coupled  in  series  to  eliminate  any 
antisymmetric  components  of  the  transient  from  the  records.  A  J" 
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diameter  steel  ball  was  dropped  from  a  height  of  4|  feet  onto  the  end 
face  of  tlus  bar  and  the  resulting  pulse  as  sensed  by  both  transducers 
was  recorded  on  an  oscilloscope.  As  shown  in  Figure  3,  certain  dis¬ 
tortions  in  the  signal  were  produced  by  the  crystal,  but  not  by  the 
surface  strain  gage,  indicating  the  need  for  corrective  Measures. 

Zn  the  next  experiment,  one  segment  of  a  similar  bar  was  drilled 

to  a  depth  of  2"  and  a  diameter  of  3/16"  as  shown  in  Figure  4  with  the 

sensor  attached  to  the  end  of  the  core.  The  core  was  then  reinserted 

into  the  hole  in  the  bar  and  the  remaining  void  was  filled  with  an 

epoxy  cement;  the  second  segment  was  attached  as  before.  A  similar 

impact  test  again  inddicated  distortion  of  the  crystal  signal;  this 

effect  was  now  believed  to  be  due  to  either  a  cross  effect  between 

longitudinal  and  lateral  strains  of  the  completely  embedded  crystal  or 

(11) 

due  to  inadequate  low  frequency  response  of  this  sensor  .  The 
latter  is  a  function  of  the  circuit  time  constant  RC  with  R  and  C  as 
the  total  resistance  and  capacitance  of  the  circuit,  respectively.  It 
had  been  estimated^)  that  if  a  3000  ppF  capacitance  were  to  be  intro¬ 
duced  in  the  circuit  by  the  addition  of  parallel  capacitors,  it  would 
permit  the  reproduction  of  both  peak  amplitude  and  duration  of  the 
pulse  to  within  at  least  95  percent.  However,  such  a  step  did  not 
materially  improve  the  signal  response.  In  consequence,  this  aberra¬ 
tion  was  thus  attributed  solely  to  the  constraint  condition  on  the 
embedded  crystal. 

In  view  of  this  conclusion,  some  attempts  were  made  to  Isolate 
lateral  crystal  surface  from  the  rock  wails.  This  was  successfully 
accomplished  by  surrounding  the  crystal  with  a  rubber  ring  that  left 


the  sides  of  the  sensor  essentially  unconstrained,  as  shewn  in  Figure 

5,  and  embedment  of  this  unit  within  the  bar  as  before.  Hie  results 

of  a  slailar  inpact  test  as  well  as  that  ox  the  previous  crystal 

arrangement,  presented  in  Figure  6  show  that  the  renoval  of  the  lateral 

constraint  on  the  crystal  pemits  the  reproduction  of  the  strain  gage 

signal  with  Mich  greater  fidelity. 

These  results  were  compared  to  those  obtained  when  charge  ampli- 

ficatlon  using  an  Endevco  Corporation  dual  channel  charge  amplifier 

with  Microdot  low  noise  cables  was  employed  to  record  the  crystal 

response  instead  of  voltage  amplification.  A  tentative  value  of  10.0 
-12 

X  10  couloobs/lb  was  chosen  for  the  piezoelectric  constant  of  the 

crystal.  Figure  7  shows  the  charge  output  for  the  constrained  and 

unconstrained  crystal  configurations,  the  latter  exhibiting  a  higher 

amplitude,  obtained  from  the  chopped  mode  of  a  single  beam  oscilloscope 

trace:  these  data  support  the  conclusiot  above  concerning  the  origin 

(12) 

of  the  distortion  in  the  signal.  As  has  been  predicted  previously  , 
the  discontinuities  due  to  bonding  of  the  two  segments  and  the  trans¬ 
ducer  embedment  had  no  observable  effect  on  the  transmission  of  the 
pulse. 

Initially,  it  was  deemed  necessary  to  check  the  consistency  of  the 

output  from  the  various  crystals  obtained  from  the  same  source  (X-tron 

Electronics,  Inc.,  Hayward);  this  was  accomplished  by  performing  a  split 

(12) (13) 

Hopkinson-bar  calibration  technique  as  shown  in  Figure  8.  X-cut 

quartz  crystals,  1/8"  in  diameter  and  1/32"  thick  were  sandwiched 
between  two  1/8"  diameter  2024-T4  aluminum  bars,  each  15"  long,  and 
subjected  to  longitudinal  impact  by  a  3/16"  diameter  steel  ball  fired 


from  a  gun  at  a  lav  Initial  velocity.  The  response  fro*,  cnis  trans¬ 
ducer  and  froa  a  pair  of  strain  gage  stations,  each  consisting  of  two 
sensors  at  the  opposite  end  of  a  dianeter  and  connected  in  series, 
located  just  ahead  and  just  behind  the  crystal  were  recorded  on  oscillo¬ 
scopes  that  were  triggered  froa  the  output  of  a  strain  gage  mounted 
about  3  inches  froa  the  impact  end.  A  Kistler  Instrument  Corporation 
Model  S/N  477  single  channel  charge  amplifier  was  used  in  conjunction 
'  with  the  crystal  Excellent  correspondence  between  the  rise  times 
and  pulse  lengths  cf  the  two  types  of  records  were  obtained.  The 
same  result  was  obtained  for  each  of  approximately  10  different  crystals. 
Indicating  tu  addition  to  the  reproducibility  feature  that  the  dis¬ 
continuity  introduced  by  the  crystal  yielded  no  measurable  effects. 

In  order  to  retain  the  lateral  freedom  of  the  crystal  upon 
insertion  in  the  specimen  by  preventing  the  bonding  material  from 
surrounding  the  sensor,  and  retain  its  freedom  from  lateral  constraint, 
it  was  found  to  be  necessary  to  mount  the  crystal  inside  a  housing. 
Initially,  this  housing  was  fabricated  from  lucite  in  view  of  its  easy 
machinability ,  and  the  unit  was  assembled  as  shown  in  Figure  9.  This 
was  then  tested  in  a  0/ 16"  split  aluminum  Hopkinson  bar  similar  to  the 
one  described  above.  The  signals  from  the  unit  showed  distinct  evidence 
of  reflections  caused  by  Impedance  mismatch  of  the  lucite  relative  to 
the  other  materials  employed.  As  a  consequence,  the  revised  design 
exhibited  in  Figure  10  was  adopted. 

Here,  the  crystal  was  encased  between  two  aluminum  end  pieces 
separated  by  a  glass  washer.  The  latter  was  cut  with  a  precision 
diamond  saw  to  achieve  parallel  faces  from  a  tube  previously  etched  to 
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the  required  dimensions.  With  this  arrangement,  the  'mpedance  mismatch 
between  the  aluminum,  quarts  and  glass  components  of  the  unit  Is 
minimised.  This  package  was  now  tested  in  the  s>lit  Hopkinson  bar 
arrangement  as  shown  in  Figure  11  with  good  results  as  portrayed  in 
Figure  12.  A  total  of  six  crystal  packages  were  assembled  and  tested 
in  this  manner;  their  response  was  in  excelle'.t  agreement  with  the 
corresponding  strain  gage  signals.  Thus,  this  technique  was  also 
employed  for  the  calibration  of  the  packages. 

A  piezoelectric  crystal  Is  essentially  a  charge  generator  and  the 
potential  difference  V  between  its  faces  is  given  by 

V  =  — -  (1) 

Ct 

where  q  is  the  charge  generated  and  ct  is  the  capacitance  of  the 
crystal.  If  the  external  capacitance  of  the  circuit  is  Cg  which 
includes  the  capacitances  of  connecting  wires,  connectors  and  the  input 
capacitance  of  the  recording  device,  then 


V  = 


c4  +  C 
t  s 


(2) 


In  the  case  of  the  X-cut  quartz  crystal  used  in  the  longitudinal  com¬ 
pression  mode,  the  piezoelectric  relations  reduce  to, 

q  =  o*  A  •  k  (3) 

where  a  i^  the  longitudinal  stress,  A  Is  the  electrode  area  of  the 
crystal  and  k  is  the  appropriate  piezoelectric  constant.  Then  Eqs. 
(2)  and  (3)  combine  to  give 


V(C  +C  ) 
t  s 

AK 


a 


vc 
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In  order  to  establish  the  values  of  the  circuit  constants  for 
the  arrangement  employed,  the  capacitances  of  the  crystal,  connecting 
wires,  connectors  and  other  external  elements  were  Measured  by  means 
of  an  impedance  bridge.  The  total  capacitance  of  the  circuit  then 
consisted  of  these  values  plus  that  of  the  plug-in  unit.  The  circuit 
resistance  consisted  essentially  only  of  the  oscilloscope  input 
impedance,  amounting  to  10  Megaohms  for  the  type  2A61  plug-in  unit 
used  in  conjunction  with  the  Tektronix  565  double  beam  oscilloscope. 
These  data  are  required  in  order  to  compare  the  piezoelectric  con¬ 
stant  determined  from  the  output  of  the  transducer  package  equated 
to  the  strain  gage  signal  with  published  values  for  this  crystal. 

The  identical  split  Hopkinson  bar  technique  used  in  the  quality 
control  studies  previously  cited  was  employed  for  the  actual  calibra¬ 
tion.  The  tests  involved  the  dropping  of  a  J  inch  diameter  sceel 
sphere  onto  the  vertical  bar  arrangement.  The  strain  gages  exhibited 
a  gage  factor  of  2.01,  a  resistance  of  120  ohms,  and  were  calibrated 

by  the  dynamic  shunting  of  external  resistances  into  the  circuit  and 

(14) 

observation  of  the  deflection  of  the  trace  .  This  permitted  the 
establishment  of  the  strain-time  record  and  the  evaluation  of  the 
transducer  package  response  upon  comparison  of  the  peak  strains,  as 
s..own  in  Figure  12.  Values  of  the  corresponding  uniaxial  stress  were 
obtained  by  means  of  t,  value  of  .0.6  x  106  psi  for  the  dynamic  modulus 
of  the  aluminum  calibration  bars. 

With  the  aid  of  Eq.  (4),  the  measured  total  circuit  capacitance, 
the  area  of  the  1/8"  diameter  crystal  face,  the  peak  voltage  output  and 
the  maximum  stress  evaluated,  the  piezoelectric  constant  was  computed 
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to  have  an  average  value  of  10.2  X  10  ±  0,2  X  10  coulombs/lb 

-12 

which  compares  well  with  the  value  of  10,01  X  10  coulombs/lb  cited 

by  Lion^^\  It  was  observed  that  the  pulses  recorded  by  the  strain 

gcges  ahead  and  behind  the  crystal  were  identical.  The  calibration 

was  performed  for  three  values  of  external  capacitance,  i.e,  3,000, 

6,000  and  12,000  upF,  and  a  total  of  4  crystal  transducers  were  tested. 

After  calibration  of  the  transducer  packages  with  a  known  external 

capacitance  utilized  to  improve  low  frequency  response,  the  v  its  were 

then  tested  with  the  aid  of  a  cha'ge  amplifier.  In  this  case,  the 

(14) 

response  is  independent  of  the  ci  cuit  capacitance  .  The  value  of 
the  piezo  lectric  constant  k  determined  from  the  calibration  tests 
was  used  to  obtain  the  charge  and  the  stress  could  then  be  directly 
computed.  The  excellent  coi -espondence  between  these  results  is 
exemplified  in  Figure  13,  A’ 1  subsequent  calibration  of  the  units  will 
utilize  the  charge  amplifier  technique  in  conjunction  with  the  estab 
1  ished  value  of  the  piezoelectric  constant. 

An  additional  physical  modification  of  the  transducer  unit  was 
required  in  order  to  permit  ready  insertion  into  relatively  long  <  ire 
holes.  This  was  accomplished  by  a  rearrangement  of  the  it  id  configura¬ 
tion  as  shown  in  Figure  14.  Although  the  acoustic  impedances  of  the 
elements  composing  the  package  are  nearly  identical,  they  differ  signi¬ 
ficantly  from  that  of  the  Yule  marble  specimen  in  which  they  are  to  be 
inserted.  However,  In  view  of  the  expectation  that  the  length  of  the 
pulses  generated  will  be  large  in  comparison  with  tue  dimensions  of  the 
transducer  package,  the  responst  is  not  believed  to  he  significantly 
affected  by  this  mismatch.  Kxperiment  s  are  curreiitli  in  progress 
utilizing  Yul1  marble  bars  to  direct lv  \crifv  this  h\po thesis. 
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Coring  of  3/16"  diameter  hole  in  Yule  marble  bars  by  means  of  a 
diamond  drill  has  shown  that  the  cores  could  not  be  recovered  due  to 
severe  crumbling  of  the  material.  Thus,  it  appeared  necessary  to  first 
drill  a  hole,  then  anchor  the  transducer  package  at  itf  base,  extract 
the  leads,  and  then  fill  and  tanp  the  remaining  void  with  a  suitable 
replacement  material.  As  a  preliminary  study,  several  1-1/16"  diameter 
bars  of  Yule  marble,  9  inches  long,  were  obtained  by  coring  out  large 
slabs  of  locally  available  material,  and  two  smooth-faced  segments 
were  cut  under  slow  hydraulic  feed.  A  3/16"  diameter  hole  about  2" 
long  was  drilled  coaxially  m  «.,ne  length  of  the  bar.  The  crystal 
transducer  W3S  inserted  and  the  hole  was  filled  with  a  mixture  of  rock 
powder  and  epoxy.  A  composite  bar  assembly  was  made  as  shown  irt  Figure 
15.  Strain  gages  ware  mounted  on  the  bar  at  the  same  section  as  the 
crystal  station.  A  dynamic  experiment  was  conduc  .ed  by  dropping  a  $" 
diameter  ball  on  one  end  of  the  vertically  held  assembly.  The  re¬ 
sulting  pulse  was  recorded  by  the  embedded  crystal  tra^lucer  and  the 
strain  gages.  Reflections,  due  to  the  epoxy  rock  mixture  filler 
material,  occurred  which  were  observed  on  the  oscillograms. 

At  present  further  experiments  are  in  orogress  to  select  a  suitable 
filler  material  which  would  closely  match  the  acoustic  impedance  of  Yule 
marble.  It  is  thought  that  this  will  consist  of  a  mixture  of  alumina 
and  epoxy  in  the  appropriate  proportions,  the  former  having  high 
impedance  and  the  latter  having  low  impedance  compared  to  that  cf  Yule 
marble.  Additional  work  is  also  underway  to  establish  a  suitable 
grouting  material.  However,  the  devel opment  work  ca  ried  out  so  lai 
indicates  that  the  installation  of  the  crystal  transducer  c  n  be 


iccompl ished  satisfactorily  in  deep  holes. 
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b)  Theoretical  Examination  of  the  Pulse  Propagation  in  an  Anisotropic 
Solid  Duo  to  a  Transient  Source 

The  formulation  of  the  equations  governing  the  boundary  value 
problem  of  a  generally  anisotropi"  half-space  subjected  to  a  normal 
time- dependent  point  loan  on  the  free  surface  are  straightforward. 

The  constitutive  relations  are 


°ij  =  cijkl  ekl  =  *  "ijkJt  (uk,A  +  “i.k* 


(5) 


where  u^  is  the  displacement  vector,  c  and  are  the  strei-  l  and 
strain  tensors,  respectively,  and  c is  tensor  the  elastic 

constants  that  satisfy  the  following  symmetry  relations 


ijk£ 


'JikJt  ~  Cijik  =  CkJiij 


(6) 


In  the  usual  tensor  notation,  a  comma  indicates  differentiation  wi„h 
respect  to  the  variable  (s)  following  it  and  the  usual  summation  conven¬ 
tion  on  spatial  variables  is  employed.  In  a  homogeneous  medium  obeying 
Eq.  (5),  for  small  displacements,  the  displacement  equations  of  motion 
under  no  body  forces  are  given  in  a  rectangular  Cartesian  coordinate 
system  xl ,  x2>  as 


Cijki  Uk,£j  ”  D  Ui,tt 


(7) 


where  o  is  the  mass  density  and  t  is  time.  When  the  free  surface  of 
the  half-space  is  given  by  x^  =  O,  the  boundary  condition  here  becomes 


Ci2kiUk,f(xr  °’  X3V  t}  '  "  6i2f(t)  6(X1?  X3) 


(8) 


where  6  is  the  Kronecker  delta,  f(t)  is  the  arbitrary  force  history 
and  6(x^,  x^)  is  the  two-dimensional  Dirac  delta  function;  the  radial 
condition  at  infinity  is  also  invoked,.  However,  the  actual  evaluation 


:  on 


of  the  fields  for  this  problem  will  require  extensive  numerical 
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computations  which,  for  the  general  anisotropic  case,  may  not  be 
feasible  at  the  present  time. 

In  the  problem  posed  above,  a  transversely  anisotropic  half-apace 
is  considered  with  the  axis  of  elastic  symmetry  located  in  the  free 
surface.  Hie  stress-displacement  equations  for  this  solid  are  given 
by 


a 


33 


e  C 


11 

Ul,l 

+  °12  U2,2 

+  C13  U3,3 

12 

U1 ,1 

+  C11  U2,2 

+  C13  U3,3 

13 

<U1,1 

+  U2,2>  + 

C33  U3,3 

44 

<U3,2 

+  U2,3) 

44 

(U3,1 

+  ul,3) 

<C11  '  '12*  <U1 ,2  +  U2,l> 


(9) 


where  the  standard  contracted  notation  c  with  m,  n  =  1,,,,,6  has  been 

mn 

employed  instead  of  the  four-index  representation  c^^jj*  "lie  solution 
of  this  problem  is  considerably  more  tractable  and  has  been  attacked 
by  two  different  procedures:  (1)  The  Integral  Transform  Technique,  and 
(2)  The  Finite  Element  Method.  Each  of  these  two  methods  has  unique 
advantages  and  disadvantages  in  the  evaluation  and  interpretation  of 
the  results. 

The  integral  transform  method  utilizing  a  Laplace-lourier  transform 

few  o\/ig\ 

and  a  Cagniard-de  Hoop  procedu  re '  is  well  established  and  permits 

interpretation  cf  the  wave  process  in  terms  of  surfaces  (velocity, 
slowness  and  wave  surface)  that  are  associated  with  phase  and  group 
phenomena  in  the  medium.  This  technique  allows  the  evaluation  of  any 
individual  response  quantity  at  any  point  of  the  medium  and  at  anv 
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specific  time  without  the  need  for  its  evaluation  at  any  other  space- 
time  point.  Thus,  the  procedure  is  very  efficient  if  only  a  few  values 
are  desired,  but  becomes  very  expensive  in  terms  of  machine  time  if 
the  entire  field  needs  to  be  evaluated  over  a  prolonged  period.  The 
numerical  methods  required  for  this  analysis  lead  to  errors  that  are 
both  known  and,  furthermore,  whose  size  can  be  controlled  readily.  On 
the  other  hand,  this  technique  also  becomes  intractable  when  the 
phenomena  to  be  investigated  involve  greater  complications  of  loading, 
anisotropy,  inhomogeneity  and  boundary  conditions.  In  this  respect, 
the  finite  element  method  has  distinct  advantages,  but  suffers  from 
the  drawbacks  of  the  need  for  the  evaluation  of  the  entire  field,  the 
lack  of  experience  in  the  use  of  this  technique  for  problems  of  this 
kind,  and  the  fact  that  the  discretization  procedure  is  not  as  well 
understood  or  controlled  as  for  the  integral  transform  analyjis. 

(1)  Integral  Transform  Technique. 

The  technique  to  be  employed  is  an  adaptation  of  that  utilized  by 
Kraut^6\  Ricketts^  and  Burridge^7\  For  the  transversely  isotropic 
case,  Eqs.  (7)  reduce  to 


°1  Ul,n  *  C2  Ul,22  *  C3  U1 ,33  +  <C1-C2,U2,21 


+  C4  U3,31  “  Ul,tt 


(C1  "  C2)  Ul,21  +  C2  U 2 , 1 1  +  Cl  U2 , 22  +  C3U2,33 


(10) 


+  C4  U3,32  =  U2,tt 


C4  Ul,31  +  °4  U2,32  +  C3  U3,ll  +  C3  U3,22  +  °5 U3; 33 =  U3 , t t 
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where  five  new  independent  constants  have  been  introduced  by 


=1  '  °U/p  c2  *  4<Cn  •  °12)/P  °3  ■  °44/P 


'4  ’  <c13  +  c44)/p 


°5  M  C33/P 


(ID 


The  equations  of  motion  and  thr  boundary  condition  will  be  subjected  to 
a  one-sided  Laplace  and  to  a  double  Fourier  transform  defined  by 


u*(x,  s)  s  |  e  9t  uk<x,  t)  dt  s  >  0,  real 


(12) 


and 


uk<!V  VV0”  JJ 


is  (O  X  +  C(,X,)  + 

e  11  **  J  u£(x,  s>  dXj  dx3 


(13) 


respectively,  where  x  is  the  radius  vector  from  the  coordinate  origin 
to  a  generic  point  in  the  medium  and  i  s/»l  ...  In  transform  space,  the 
former  become 

[c2D2-s2(cl<yJ+c3<y^fl)]  5J  -  (cl-c2)uori  J  ‘  c4  "  0 

"  (C1  “  C2)i8®lCGl  +  tcl°2  "  s2<c2ffl  4  c3°3  +  ‘  c4  ls0f3  D“3  =  0  (14) 

•'(‘Vs"!  -  c4iW3  *  rc3D2  -  '’<c3 a'i  *  C5°3  *  0 


where  the  symbol  D  denotes  differentiation  with  respect  to  x  .  The 
corresponding  boundary  conditions  on  =  0  are 

Du*-  isQf^u*  =  0 


<C1  ~  2c2>  iso^u*  -  CjDu’  +  (°4  -  C3)  isa.^  =  f*(s)/p 


_* 


(15) 


DU3"  iSV2  =  ° 


A  solution  is  now  assumed  in  the  form 

„  _S^X2 
\  ■  \  0 


(16) 
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which,  upon  substitution  in  Eqs.  (14),  leads  to  a  system  of  homogeneous 
algebraic  equations  in  U^;  the  condition  for  a  non-triviel  solution 
yields  the  following  deterainantal  equation 


^C2X2'Cl"l  "  C3ff3  "l) 


(cl-c2)l«iX 


-  c.a,er_ 
4  13 


<Cl”C2)lflfiX 

<ClX2_C2(>l‘C3a3‘l) 


'  Ws 


C4l0f3X 


W 


-  °3ai~  Vs'l>i 


«  0  (17) 


Equation  (17)  is  a  sextic  in  X  and  defines  the  slowness  surface.  For 
a  transversely  isotropic  material,  it  degenerates  into  a  quadratic  and 
a  quartic  factor  which  greatly  simplifies  the  subsequent  analysis. 

The  slowness  surface  is  given  by 
[(CjA2  -  II)  (c3A2  -  N)  +  c20f2  A2]  (cgA2  -  M)  =  0  (18) 

where 

M  =  c3o2  +  1 ,  N  *  CgCr2  +  1 ,  and 


Consequently,  the  expressions  for  X  *  X.(a  ,  w„),  with  j  s  1,  2,  3 

J  J  i  O 

say  be  written  as 


X1  "  *2°  =  K  +  +  l)/c2]3 


*2,3  ~  *2 


±Ji2  -  4c^c3KN 


2c,  c 


1L3 


(19) 


2  2 

where  i  =  c^N  +  c^M  -  0f3  and  X^  are  to  be  chosen  such  that  Re(X^)ao 

to  satisfy  the  radiation  condition.  Substitution  of  these  eigenvalue's 


into  the  homogeneous  algebraic  equations  yield  the  three  associated 
(j) 

eigenvectors  Uk  so  that  the  general  solution  to  Eq.  (14)  now  takes 


the  fora 


K-  l  » 


(!)  e‘‘XJX2 


It  the  three  Independent  are  choaen  as  y^\  then  Eq.  (20)  nay  be 


written  as 


;*  .  U«>  a"'*1’'2  ♦  „<»  e"^2  e  u<3>  a"’*3*2 
“J  =  -  (IOj/Xj)  Oj11'  e  1  2-  W/V  °1<2>*  '^l'2 
-  (tyaj)  U<3)  a"’’'3*2 

=5  -  sV  o,  »;2>  • 5X2X2  *  »3  »{3)  «“V2] 

X  o 


where 


bj  =  (c^  -  M)/e4  ;  j  s  2,3. 


The  general  solution,  Bq.  (21),  is  then  adjusted  to  satisfy  the 
boundary  conditions,  Eqs.  (15),  yielding  three  equations  in 

<X1  +  <y?)0l<l>  +  2\\  Ul<2>  +  2X1X3  Ul<3>  *  0 

+  L2U12^  +  L3U1^  =  (22 

Vs  Ul(1>  +  XlX2(Q'3  +  b2)Ul<2>  +  W°S*  V0^  =  ° 


where 


L1  85  2c2V 


hsll+cl^-  <c4  -  c3)bk  ’  k  =  2-  3* 

The  relations  are  now  solved  for  the  U*J),  and  substitution  of  these 
values  into  Eq.  (21)  yields  the  d ispl acement  field  in  transform  space 
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j,  *  lflri  r  “*xix2 
“i  *  '  <•>  &Tn  Lai  * 


'■Va 

+  A2  e  ■*•  Ag  e 


-•Van 


1  r“  1  “SX  X_  "sX2X2  -sX3Xo“l 

u2  - f  <8)  pi:  [r  V  +  V» 6  +  Vs 6  J 


0  1 


-s\„x. 


-sX„x, 


-*  „  1  r.  -  “'2"2  .  .  ”'3"2'\ 

=1  ■  f  <*>  wTdT  LVa  e  *  baA3  e  J 


3  0 


(23) 


where 


Dq  »  AjL^  +  AgLg  +  A3L3 


A.  »  -  2X  XX_  (b  -  b  ) 

X  1  •  3  <3  « 

A2  *  W0^  +  V  +  2Q?1  b3-^ 

A3  «  -  X2[A2(cy2  +  b2)  +  2a2  b2] 


Us it.:  the  Fourier  inversion  theorem  yields  the  solution  in  the  Laplace 
transform  space  as 


®  3 


-s ( icy  x,  +  io0x_  +  X.x„) 


■:<*.■> -“.vv  !i"'n'T2, 

4tt  p 


j=l 


d(Xlda3 


(24) 


where 


w(j)  =  ^1 


1  Dq  "j 

2 


A ,  ,  J  a  1,  2,  3 


W(1)  -  j!L  A  W(J)  -  i»  1-23 
W2  "^Al  '  W2  *  Dq  AJ  ’  J  "  2>  3 


*31>  "  0 


(j) 

•V  =^V'2-3 


The  inversion  cf  Eq.  (2-1)  into  the  real  time  domain  is  expedited 
by  the  transform 


=  w  cos  9  =  -  ip  cos  9 
Qfg  a  w  sin  9  a  -  ip  sin  9 


-  «  <  w  <  06  ;  -  Jtt  *i  9  <  Jn  (25) 
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where  the  Jacobian  of  the  transformation  ia  |j(p,  8)J  =  p  and  p  is 
considered  to  be  complex.  Furthermore,  a  polar  coordinate  system  in 
the  x^Xg-plane  is  introduced,  given  by 

x,  a  r  cos  x„  *  r  sin  with  0  r  <  «°  and  0  £  8  <  2n  (26) 

1  3 


Equation  (24)  is  now  rewritten  as 


TT 

2  i» 


-sfprcos  (K-8)+X  x_] 


uj<x,s>. - VJ de  J  •*•«•>  I  . ~J***  ») 


4tt  p 


TT  -  j[os 
"2 


!=1 


(16) 

It  is  inverted  utilizing  the  Cagniard-de  Hoop  method 
variables  are  introduced  by  means  of  the  relations 
Re[pr  cos(S-Q)  +  Xj(p,8)  x2]  =  t^  , 
«9m[pr  cos(1)-8)  +  Xj(p,8)  xg]  =  0 


in  which  new 


(28) 


where  t  is  real  and  positive.  Hiese  relations  define  three  sets  of 
J 

parametric  relations  in  t  for  the  Cagniard-de  Hoop  path  on  the  complex 

p-plane.  Using  Cauchy's  Integral  theorem  and  Jordan’s  Lemma,  it 

can  be  shown  that  integration  along  the  Bromwich  contour  for  Eq,  (27) 

(6) 

is  equivalent  to  integration  along  the  Cagniard-de  Hoop  pathv  , 
Equation  (28)  may  be  expressed  as  a  polynomial  p  as 

2f  g  2  2  c3  2  2 

p  {cos  cp  cos  (8-9)  +  (cos  0  +  —  sin  9)  sin  tp } 

C2 


2  <29) 
-  p[2T  coscp  cos  ($-8) ]  +  [T2  -  1±2-SL]  =  o 

C2 


and 

P4+ A(T,'§,‘5,cp)p3+ B(T,‘5,8,Cp)p2  +  ClT.’S.e^Jp  +  E(T,^,8,cp)  =0  (30) 

for  the  quadratic  factor  involving  X^  and  the  quartic  factor  involving 
\2  and  X^,  respectively,  in  Eq.  (19) 
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where,  T  =  t/R, 

A  =»  -  ~  (2t  cos3<p  cos3($-0)  +  TG  coe^  cob($-8)  8in2q>]  , 
F 

B  *s  -  [6t2  cos2(p  cob2(§-8)  ♦  t2G  sin2tp  -  sln2f  cos2cp  . 
cos2  $-8)  ©j  -  sin4tp  J] 

C  =  -  ||  [2t3  coscp  cobC5-8)  -  T  cos  cp  sin2®  COS$f  8)©^] 

E  =  ^  [T4  -  T2  sin2cp  ©j  +  sin4cp  h^] 

F  =  cos4<p  cos4  $-8)  +  sin2cp  cos^p  cos2(5-0)G 
+  sin4(p  K 

K  a  cos48  +  (2-m.)  sin28  cos28  +  —  sin40 
1  cl 


2  2 

J  =  cos  0  +  sin  0  , 


G  =  2  -  sin  0 

2  2 

m,  a  h,  (2c, c_  -  c,c  -  c_  +  c.) 

1  I  1  o  10  o  4 


ei  •  hi<ci  *  'a* 


£1  =  hl<C3  +  C5} 


gl  =  hl<Cl  '  C5) 


hl  =  1/C1C3 


In  the  above,  an  additional  space  transformation  is  introduced,  given  by 

R  =  (r2  +  x2)^  ,  0  s  r  <  co  ;  cp  =  tan  1  —  ,  0  £  cp  <  n  (31) 

&  r 

Since  Eqs.  (29)  and  (30)  are  polynomials  in  powers  of  p  with  real 
coef f U’ients ,  their  roots  are  either  real  or  complex  conjugate.  Upon 
deforming  the  path  of  integration  of  Eq.  (27)  into  Cagniard-de  Hoop 


contours,  recalling  that  s  and  t  are  real  and  taking  advantage  of  the 
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symmetry  of  the  path  with  respect  to  the  real  axis  of  the  p  plane, 
there  results 

TT 


2  •  £ 

u*(x,s)®~~—  Jd6  |  sf*(s)  ^  [Re  (p,e,$,tp) >  ||  e  J  dt 

2"PE  «  J=1 

2  J 


(32) 


The  inversion  of  this  relation  into  the  real  tine  domain  is  now  obtained 
hy  inspection  as 


TT 

2  t 


«k<x.t)  =-~ —  JaeJ  l(t-t')  ^  <U[w'J><p,e,e,<e)] 

*"  PR  -  *  j.i 


TT  t. 


p(t')  dt' 


(33) 


To  effect  this  procedure  with  any  degree  of  facility,  it  was  necessary 
to  define  Rienann  surfaces  for  the  Cagniard-de  Hocp  paths  defined  by 
Eqs.  (29)  and  (30)  where  the  values  of  p  are  single-valued  and  analytic. 
Since  there  are  three  branches  of  X  derived  from  the  sextic  equation 
(19),  there  must  exist  three  such  Riemann  sheets.  If  such  surfaces  are 
not  utilized  in  the  inversion  process,  then  p  will  be  a  multi-valued 

(6) 

function  and  the  analysis  becomes  vastly  more  complicated  and  tedious' 
The  time  function  f(t)  will  now  be  replaced  in  Eq.  (33)  by  the 
Dirac  delta  function  6(t),  and  the  Heaviside  function  H(t)  =  *<q  * 
respectively,  yielding 

TT 

2  3 

V;,t>='  rh;  ^  d8[s?{X  ||  wt-t  >]  (34) 

0  2TT  pH  j  t  ^ 
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n 

2  3 

u  =  -  -i—  J  d«[  2,«et»*J)(p.«,?,ip)J  p  |S  Wt-t  )j  (35) 

H  ""“.p  J,l 

"  2 

Equation  (34)  is  the  Green's  function  for  the  present  dynamic  problem 
and  the  displacement  field  for  an  arbitrary  time -dependent  force  f(t) 
may  be  obtained  by  convolution  of  this  function  with  f(t).  Equation 
(35)  represents  another  relation  that  may  be  employed  to  obtain  a 
general  solution  for  an  arbitrary  time  history  upon  convolution  with 
~  f(t).  The  latter  is  a  more  convenient  relation  for  this  purpose 
since  the  kernel  function  need  not  be  differentiated  with  respect  to 
time. 

It  should  be  noted  that,  although!?  and  6  are  independent  variables 

in  two  separate  spaces,  the  integrands  of  Eqs.  (33)  to  (35)  are  related 

to  the  response  on  the  =  constant  plane  due  to  a  corresp  1  in.g  line 

source  at  0  =  0  90  in  the  two-dimensional  solution  due  to  the 

principle  of  duality  between  the  slowness  space,  where  8  is  defined,  and 
(17) 

the  wave  surface  .  Consequently,  the  solutions  of  these  equations 
can  also  be  considered  as  the  superposition  of  the  responses  of  a  con- 
ti  <?  distribution  of  line  sources  in  all  directions  within  the  x^  =  0 
plane  upon  multiplication  by  a  weighting  factor. 

For  the  propagation  of  surface  waves,  the  Cagniard-de  Hoop  path 
collapses  to  the  real  axis  of  p  which  contains  the  huyleigh  pole.  Since 
the  path  of  integration  must  be  deformed  to  exclude  this  singularity, 
the  resulting  infinitesimal  semicircle  no  longer  corresponds  to  real 
time  and  the  contribution  from  the  Rayleigh  pole  must  be  evaluated 
separately.  This  contribution  is  given  by 


,  y^.»«33S |rtf 
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<x,,  0,  x,.  t)  =  R  •  6C  t  - - - )  (36) 

“pole  1  3  1  yS,0) 

where  v  (9,0)  is  the  Rayleigh  wave  velocity  obtained  f row  the  location 

R 

of  the  pole  and  constant  is  given  by 

TT 

2 

B1  =«e  I  J  Q(9,  6,  p)  d6|  (3.) 

n 

"  2 

with  Q  as  the  residue  of  the  expression 

3 

£  M  ^j>  (p,  6,  *51  (38) 

J=1 

appearing  in  Eqs.  (34)  and  (35)  at  the  Rayleigh  pole. 

Except  for  the  deterainat ion  of  the  stress  field,  the  equations 
presented  above  conclude  the  formal  analysis  of  the  problem  and  the 
evaluation  of  the  field  variables  requires  the  use  of  appropriate 
numerical  methods  as  detailed  below.  It  is  expedient  as  an  intermediate 
step  to  compute  certain  physical  features  closely  associated  with  the 
analysis,  i.e.  the  slowness  and  wave  curves  and  the  Cagniard-de  Hoop 
paths  for  the  body  waves. 

Yhe  slowness  curves  are  evaluated  by  fixing  either  6  or  9  in  Eq . 

(17),  Figure  16  shows  a  set  of  three  slowness  curves  for  a  value  of 

9  =  99*  with  varying  from  1°  to  89°  in  increments  of  2°,  with  the 

s^  and  axes  corresponding  to  the  inverse  of  the  phase  velocities 

along  the  x^  and  x^  axes  of  the  medium.  Tie  dynamic  elastic  constants 

( 8 ) 

employed  in  this  analysis  are  those  utilized  by  Ricketts'  as  shown  in 


Table  1. 
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TABLE  1.  VALUES  OF  DYNAMIC  EU  STIC  CONSTANTS 

FOR  A  YULE  MARBLE  BLOCK 

10  2  2 

c1  =  3.8  X  10  in  /sec 

10  2  2 

c2  =  1.28  X  10  in  /sec 

10  2  2 

c  =  0.745  X  10  in  /sec 

*3 

10  2  2 

c.  =  1,89  X  10  in  /sec 
4 

10  2  2 

c  =  3,33  X  10  in  /sec 
d 

Such  curves  have  already  been  reported  in  Reference  (8)  for  both  body 
and  surface  waves. 

Wave  curves  can  also  be  obtained  by  fixing  either  6  or  X  and  are 
usually  evaluated  by  the  rather  complicated  procedure  o.  taking  polar 
reciprocals  of  the  slowness  surface.  In  the  present  analysis,  they 
will  be  determined  by  a  simpler  procedure.  By  virtue  of  Kq.(28),  the 
Cagniarl-de  Hoop  transformation  may  be  considered  as  the  equation  of  a 
plane  in  the  slowness  space  which  is  normal  to  the  position  vector  of 
receiver  station.  With  increasing  time,  this  plane  moves  outward 
from  the  origin  while  maintaining  its  normality  with  the  position 
vector.  At  some  particular  time  t^  ,  this  plane  will  become  tangent  to 
a  slowness  surface  and  it  can  be  then  considered  as  a  point  on  the  wave 
surface  from  the  principle  of  duality.  Thus  the  wave  surface  an  be 
regarded  as  being  composed  of  the  totality  of  such  tangent  planes,  each 
associated  with  a  different  radius  vect.  r. 

Furthermore,  these  tangent  points  correspond  t.o  the  double  roots 
of  the  t’agn  l  a  i’ll -do  Hoop  paths  defined  b>  tqs.  (29)  and  (3<>).  Thus  the 
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discriminants  D  associated  with  the  solutions  for  p  of  these  equations 
P 

must  vanish  and  the  roots  are  real  since  the  coefficients  in  the 
equations  are  real.  Consequently,  the  first  of  these  conditions  permit! 
the  evaluation  of  the  quantity  t,  the  arrival  time  at  a  given  receiver 
station,  from  a  relation  of  the  form 

Dp(T,  0,  6,  cp)  =  0  (39) 

By  M.rodueing  the  varifble  r  =  1/v  l».  Eq.  (39)  and  solving  this 
equation  for  various  values  of  p,  each  corresponding  to  a  value 
v^0 ,  0 ,  cp)  ,  a  plot  of 

Ri  =  Vi >  6,  cp)  *  t  at  t  =  1  (40) 

provides  the  wave  surface.  Figures  17  and  18  show  the  pair  of  wave 

curves  on  the  plane  It  -  90  for  cp  varying  from  1°  to  90°  and  on  the  plane 

<p  =  45°  for  various  values  of  which  are  obtained  from  the 

Dp(T,  0,  8,  cp)  =  0  of  Eq.  (30),  respectively.  The  first  of  these 

diagrams  indicates  that  there  exist  three  characteristic  ranges  for 

this  medium,  0  -s  cp  <"  cp^,  cp  cpp  <  rt>2  and  cp2  *-  cp^  <  90°,  while  the 

second  portrays  two  ranges  for  0,  namely  0  <  6  <  $  and  ^  <-  0p  '  90°, 

The  cusps  in  these  curves  represent  characteristics  not  exhibited  by  an 

isotropic  medium,  Thu3  the  normal  number  of  roots  of  Eq.  (39)  is  three*, 

but  in  the  critical  angular  ranges  Cp^  and  Tf  ,  there  are  four  such  roots 

for  a  transversely  isotropic  material, 

A  computer  program  was  developed  which  solves  the  quart ic  equation 

(18) 

(30)  denoting  the  Cagniard-de  Hoop  path,  employing  Kairstow's  method' 
and  the  results  are  indicated  in  Figures  19,  20,  and  21.  The  curves 
clearly  indicate  the  nature  of  the  critical  regions  of  the  wave  surfaces 
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T? 

4 


as  discussed  above,  with  only  Figure  20  showing  four  intersections  of 
the  path  with  the  real  axis.  Sheets  I,  II,  and  III  of  the  path 
correspond  to  values  of  X^  ,  X^  and  X^,  respectively,  each  converging  to 
an  asymptote  for  large  values  of  t.  A  number  of  other  Cagniard-de  Hoop 
paths  and  slowness  and  wave  curves  have  been  evaluated,  but  are  not 
reproduced  here  since  they  do  not  exhibit  any  additional  unusual 
features, 

A  computer  program  has  also  been  written  to  evaluate  the  displace¬ 
ment  field  from  Eq.  (35);  this  is  currently  in  the  process  of  debugging. 
It  may  be  necessary  to  revise  and  improve  the  program  in  order  to 
reduce  the  required  computer  time  for  the  location  of  the  real  and 
complex  roots  of  the  Cagniard-de  Hoop  path,  a  process  that  appears  to 
be  more  costly  than  necessary.  The  program  also  determines  the  real 
roots  for  Eq.  (39)  by  equating  the  discriminant  to  zero  and  entails  a 
polynomial  evaluator  and  an  integration  scheme  employing  Simpson's 
rule.  Presently,  a  period  of  two  seconds  is  required  to  compute  the 
displacement  field  at  a  specified  point  and  a  given  instant  of  time. 

The  total  computer  time  for  the  evaluation  of  a  displacement  fieid 
along  a  given  ray  of  the  domain  resulting  from  a  Heaviside  input  is  a 
linear  function  of  the  number  of  time1  intervals  that  are  to  be  included 
in  history  of  the  transient.  It  should  be  noted  that,  for  a  given  rnv , 
the  displacements  are  calcuated  for  given  values  of  i  It,  and  thus  the 
history  at  all  points  along  such  a  line  can  be  scaled  m  accordance  with 
this  parameter.  For  fifty  time  intervals,  approximately  one-half  hour 
of  computer  time  on  the  CI)C  6600  would  be  required  to  evaluate  the 
displacements  along  50  rays.  In  addition,  a  convolution  integral 
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procedure  must  be  programmed  in  order  to  accommodate  specified  input 
histories;  such  a  process  is  estimated  not  to  be  very  time-consuming. 

Figure  22  presents  preliminary  results  for  the  Cartesian  dis¬ 
placement  components  u1 ,  Ug,  u3(x,  t)  for  ^  =  80°  and  9  =  45°  for  a 

-5 

sequence  of  19  time  intervals,  beginning  with  a  spacing  of  0.1  X  10 

-5  -5 

sec/'n  from  t  =  0.4  X  10  to  0.8  X  10  sec/in  and  an  interval  of 
-5 

0.2  X  10  sec/in  beyond  this  range.  The  time  of  computation  was  about 
J7  seconds;  the  angle  9  =  45°  was  chosen  since  it  lies  in  the  critical 
range.  This  diagram  exhibits  drastically  different  characteristics 
than  for  the  corresponding  two-dimensional  solution'  ,  primarily  in 
the  absence  of  strong  discontinuities  that  are  present  in  the  two- 
dimensional  case.  Displacements  in  cylindrical  and  spherical  components 
have  also  been  computed  for  subsequent  comparison,  but  are  not  included 
in  the  present  report, 

(2)  Finite  Element  Method. 

A  finite  element  computer  code,  FEAP-71  developed  by  Professor  R. 

L.  Taylor  of  the  Department  of  Civil  Engineering,  University  of 
California,  Berkeley,  is  being  adapted  for  use  on  the  subject  problem. 

It  is  a  research  oriented  finite  element  assembly  program  with  a  current 
selection  library  of  twenty  elements  and  extensive  input-output  utility 
routines  and  automatic  error  checking.  Extensions  developed  here  con¬ 
sist  of  the  construction  of  a  three-dimensional  element  and  associated 
time  integration  schemes  suitable  for  use  in  a  wave  propagation  problem. 

The  dynamic  formulation  of  the  finite  element  method  can  be  derived 
from  Lagrange's  Central  Principle  which  combines  virtual  work  and 
D'Alembert's  Principle.  The  virtual  work  statement  in  a  continuous  body 
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can  be  expressed  in  indicial  notation  as 

J  6ui  dV  +  i  J  aii  6eii  dV  “  I  Tt  6ui  dS  +  J  Fi  6ui 


dV 


(41) 


where  6  is  the  variation,  and  T^,  V  and  S  are  the  components  of  the 
surface  traction  and  body  force,  and  the  volume  and  surface  of  the 
body,  respectively.  In  order  to  obtain  the  displacement  equations  of 
motion,  the  stress  and  strain  terms  in  Eq.  (41)  are  converted  by  means 
of  the  constitutive  equation  and  the  strain-displacement  relations; 
the  surface  integral  is  transformed  into  a  volume  integral  by  means  of 
the  divergence  theorem.  The  finite  element  approximation  is  now  intro¬ 
duced  by  subdividing  the  domain  V  and  introducing  a  local  Ritz  technique 
such  that  for  each  element,  the  displacements  are  approximated  by 

u ^ ( x  ,  t)  =  (  tp^x)  )  fu(t)1  (42) 

where  (cp^x))  is  the  row  vector  of  prescribed  interpolation  functions 
and  fii(t)}  are  the  nodal  displacements  represented  by  a  column  vector 
of  generalized  coordinates.  Insertion  of  Eq.  (42)  into  the  displace¬ 
ment  version  of  the  virtual  work  statement  yields  the  set  of  coupled 
1 tnear  differential  equations 


[M]  {u}  +  [K]  {u}  a  [f] 


(43) 


where  [M]  is  the  diagonal  mass  matrix,  [K]  is  the  global  stiffness 

matrix,  [u]  is  the  nodal  displacement  vector  fo^  the  encire  domain,  {ul 

is  the  nodal  acceleration  vector,  and  {fl  is  the  loading  vector. 

The  solution  of  Eq.  (43)  is  accomplished  by  using  Newmark's  Beta 
(19) 

Method  which  is  a  step-by-step  forward  integration  procedure,  given 
by 


l-5  r" 
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u  ,  a  u  +  (1  -  V)  u  At  +  Yu  .At  (44) 

n+1  n  n  n+1 

u  ,  =  u  +  u  At  +  (£  -  8)ti  (At)2  +  0u  .(At)2  (45) 

n+1  n  n  n  n+i 

where  8  is  a  parameter  characterizing  acceleration  over  the  interval 
At,  V  is  a  quantity  that  introduces  artificial  damping,  At  is  the  time 
step  increment  and  ur,  u;  and  are  the  displacements,  velocities 
and  accelerations  at  time  step  n,  respectively. 

Upon  inspection  of  Eqs,  (s3)  -  (45),  it  is  clear  that  the  solution 
algorithm  is  implicit  in  that  the  displacement  at  time  t  +  At  is 
dependent  upon  the  acceleration  at  this  time,  thus  requiring  either  a 
costly  direct  solution  or  an  iterative  procedure.  However,  if  the 
Newmark  parameter  8  is  chosen  as  zero,  an  explicit  algorithm  results 
so  that  the  displacement  at  the  new  time  can  be  obtained  from  infonia- 
tlon  available  at  the  previous  instant.  Such  a  choice  requires  the 
exercise  of  considerable  care  to  insure  that  the  numerical  procedure 
is  stable.  Experience  has  shown  that  stability  can  be  achieved  by 
appropriate  selection  of  the  time  step,  namely  that  the  time  step  be 
less  than  the  ratio  of  the  minimum  inter-node  distance  to  the  maximum 
wave  speed.  This  approach  was  selected  here  to  substantially  reduce 
computer  time. 

The  majority  of  the  effort  in  this  phase  thus  far  has  been  the 
coupling  of  the  Newmark  explicit  method  with  the  original  finite  element 
code,  FEAP-71,  In  addition,  several  programming  tricks  have  been 
incorporated  which  greatly  reduce  computer  storage  and  time  requirements. 
The  half-3pace  is  modelled  by  a  mesh  of  identical  element  i  that 
requires  the  formation  of  only  a  single  element  stiffness.  Moreover, 
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the  formation  of  the  global  stiffness  matrix  [K]  which  requires 
massive  storage  on  slow  speed  devices  has  been  abandoned.  Instead, 
the  operations  involving  this  matrix  from  Eq.  (43)  are  performed  by 
manipulating  single  element  stiffness. 

After  assembly  of  this  program,  initial  efforts  were  concentrated 
on  the  debugging  phase  in  which  a  simple  test  problem  was  devised  that 
consisted  of  the  axial  loading  by  a  triangular  pulse  of  a  stack  of 
cubic  elements  with  eight  nodes  per  element,  termed  Brick-8.  This 
three-dimensional  model  collapses  to  the  uni-axial  homogeneous  elastic 
wave  propagation  problem  when  the  medium  is  treated  as  a  homogeneous, 
isotropic  elastic  substance  with  a  Poisson  ratio  of  zero.  The  results 
of  this  check  are  presented  in  Figure  23  which  shows  the  excellent 
correspondence  of  the  predict  ions  of  the  finite  element  method  and 
that  of  classical  analysis. 

The  confidence  gained  by  this  agreement  prompted  the  construction 
of  a  more  complicated  test  problem,  close  in  character  to  the  case 
under  consideration  in  this  study,  for  which  closed-form  analytical 
solutions  as  well  as  other  numerical  data  were  available.  This  consisted 
of  the  determination  of  the  response  of  the  isotropic  elastic  half-space 
under  concentrated  time-dependent  normal  load,  as  shown  in  Figure  24a, 
with  the  three-dimensional  finite  element  model  shown  in  Figure  24b. 

Here,  the  most  difficult  aspects  of  the  entire  procedure,  namely  the 
choice  of  a  suitable  element  including  higher-order  interpolation  func¬ 
tions,  the  time  step  and  the  spatial  discretization  were  to  be  examined 
initially.  The  region  in  which  it  is  most  difficult  to  obtain  accurate 
results  by  the  finite  element  method  is  the  domain  close  to  the 
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concentrated  load  and  on  the  surface  of  the  half-space  where  dis¬ 
continuities  associated  with  the  Rayleigh  wave  dominate.  If  suitable 
correlation  of  the  finite  element  computations  with  closed-form 
solutions  are  obtained  in  this  region,  then  even  better  agreement  is 
expected  to  exist  in  the  remainder  of  the  space. 

In  view  of  the  fact  that  extensive  testing  of  the  efficacy  of  the 
clement  and  optimality  of  the  spatial  discretization  on  a  three- 
dimensional  basis  is  extremely  expensive,  this  aspect  of  the  numerical 
investigation  was  performed  initially  using  an  axisymmetric  (i.e.  two- 
dimensional)  finite-element  formulation  as  indicated  in  Figure  24c. 

'Ibis  code  is  adequate  to  ascertain  the  validity  of  the  choices 
Involving  these  two  items,  but  it  cannot  be  employed  in  the  analysis 
of  the  actual  problem  since  the  anisotropy  encountered  here  destroys 
the  axisymmetric  nature  of  the  phenomenon.  Thus,  upon  completion  of 
this  phase  of  the  numerical  investigation,  the  three-dimensional  code 
will  be  tested  directly  with  particular  emphasis  on  the  aspects 
relating  to  the  anisotropic  nanre  of  the  problem. 

The  initial  results  involving  a  comparison  of  the  finite-element 

method  and  a  numerical  convolution  with  a  closed-form  solution  of  this 

(20)  2 
problem  for  a  pulse  of  the  form  f(t)  --  (sin  2rrt)  for  0  ‘  t  *'  0.5 

and  zero  otherwise  applied  at  the  coordinate  origin  to  the  surface  of 

an  isotropic  half-space  with  a  Young's  modulus  K  --  2.5,  a  Poisson  ratio 

v  =  $  and  a  density  p  •-  1  are  presented  in  Figures  25a,  25b,  and  25c. 

The  shape  of  the  pul.se  corresponds  to  that  observed  experimentally  at 

( H ) 

the  impact  point  when  a  steel  sphere  is  dropped  on  a  block  of  rock 
The  diagrams  exhibit  the  vertical  surface  displacements  at  three 
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positions  near  the  impact  point  as  a  function  of  time  for  both  the 
three-dimensional  and  the  axi-symmetric  finite  element  analysis,  as 
well  as  the  results  derived  from  the  closed-form  solution.  It  may  be 
observed  that,  although  the  correspondence  between  the  three  sets  of 
data  nearest  the  contact  point  is  not  outstanding,  it  dramatically 
improves,  as  expected,  with  distance  from  the  impact  point  as  would 
also  be  anticipated  to  occur  in  the  interior  of  the  medium.  In  parti¬ 
cular,  there  is  close  correspondence  between  the  two-  and  three- 
dimensional  finite  element  results  at  points  x  =  0.2  and  0.3  inches, 
indicating  the  validity  of  the  previous  conclusion  that  extensive 
testing  of  the  efficacy  of  the  code  can  be  performed  with  the  simpler 
model.  As  further  modifications  of  the  program  are  initiated,  it  is 
expected  that  the  concurrence  between  the  data  presented  here  will 
improve. 

c)  Property  Determination  of  Yule  Marble 

The  overall  objective  of  this  phase  of  the  investigation  consists 
of  the  determination  of  the  geometric  and  mechanical  properties  of 
Yule  marble,  including  fracture  properties,  as  a  function  of  strain 
rate.  The  material  is  considered  to  be  transversely  isotropic  based 
on  both  previously  published  information  concerning  this  rock  and  the 
nature  of  its  geologic  formation,  and  its  crystal  size  is  such  that  it 
can  be  considered  as  macroscopica' ly  homogeneous  in  relation  to  the 
wave  lengths  of  transients  produced  in  presently  available  experimental 
specimens.  Five  major  subdivisions  of  the  testing  program  can  be 
differentiated:  (1)  The  determination  of  the  axis  of  transverse 

isotropy  of  the  material,  principally  by  use  of  crystallographic 


techniques,  (2)  the  execution  of  quasi-static  tests  on  specimens  whose 

axes  are  parallel  or  perpendicular  to  the  axis  of  transverse  isotropy 

to  provide  stress-strain  curves  in  tension  and  compression,  at  rates 
-5  -1 

ranging  fror„  10  to  10  per  second,  as  can  be  obtained  in  standard 
testing  machines,  (3)  the  development  of  procedures  to  obtain  inter¬ 
mediate  strain  rates  of  the  order  of  10  to  100  per  second  by  modifying 
the  sp,.t  Hopkinson-bar  technique,  (4)  the  specification  of  the 
variation  of  the  behavior  of  the  materiel  in  tension  and  compression 
as  a  function  of  strain  rate,  and  (5)  the  characterization  of  the 
failure  process  in  terms  of  the  strain-rate  parameter,  including 
optical  observations  of  the  phenomenon.  From  this  information,  a 
comprehensive  model  of  the  mechanical  response  of  the  substance  may 
be  evolved. 

The  determination  of  the  axis  of  transverse  isotropy  is  accom¬ 
plished  crystallographically  by  microscopic  examination  of  thin  sections 
obtained  from  three  mutually  orthogonal  direct  ions.  These  are  prepared 
by  first  cutting  the  rock  into  thin  slabs,  attaching  these  slices  to 
glass  slides  with  Canada  balsam,  and  then  lapping  the  rock  with  fine 
abrasive  until  the  desired  thickness  is  obtained.  Since  marble  is  a 
matrix  of  calcite  crystals  with  optical  axes,  oriented  randomly  about  a 
general  preferred  direction,  it  is  possible  to  determine  each  part icular 

crystal  orien  ntion  by  measuring  the  angle  between  the  optic  axis  of  a  given 

(21 ) 

crystal  and  some  fixc-d  reference  direction  and  plotting  this  rela¬ 
tionship  on  an  equal  area  projection.  When  a  large  number  of  crystal 
axes  are  located,  a  general  pattern  of  average  crystal  orientation 


erne rges. 
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These  measurenents  are  performed  by  locating  a  uniaxial  inter- 

(21) 

ference  figure,  a  metalope  ,  produced  by  light  passing  along  the 
optic  axis  of  a  calcite  crystal.  Each  crystal  is  observed  under  a 
microscope  equipped  with  a  universal  stage  which  permits  the  specimen 
to  be  rotated  about  two  independent  axes  corresponding  to  Euler 
angles  so  that  the  optic  axis  can  be  aligned  to  produce  a 
metalope  centered  on  the  crosshair. 

While  X-ray  techniques  provide  another  way  of  establishing  the 
orientation  of  the  axis  of  transverse  isotropy,  this  method  is  elab¬ 
orate,  time-consuming  and  expensive  including  the  necessity  of  some 
computer  programming.  In  consequence,  it  is  not  anticipated  that 
this  procedure  will  require  employment  in  the  present  investigation. 

In  order  to  retain  the  integrity  of  the  two  large  marble  blocks 
intended  for  wave  propagation  studies,  slabs  with  dimensions  of  1'  X 
2j'  X  lj"  were  purchased  from  Clervi  Marble  Company,  San  Francisco 
for  the  purpose  of  furnishing  specimens  for  static  testing.  Since  the 
axis  of  transverse  isotropy  for  this  slab  was  unknown,  cylindrical 
test  samples  3/4"  in  diameter  X  2"  long  were  simply  cored  from  this 
source  in  three  mutually  orthogonal,  albeit  arbitrary  directions  and 
lapped  to  within  5/10, 000th  of  parallel.  In  order  to  prevent  crumbling 
or  splitting  out  during  the  coring  process,  a  glass  plate  was  cemented 
with  Canada  balsam  to  the  distal  side  of  the  slab. 

The  state  of  strain  in  these  specimens  were  measured  by  means  of 
SR-4  foil  resistance  gages  mounted  by  means  of  Epoxy  150  cement.  A 
compromise  was  required  to  select  a  suitable  size  for  this  transducer. 

A  gage  of  excessive  size  records  the  avenge  strain  across  a  region,  not 
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at  a  point,  and,  further,  is  more  difficult  to  attach;  a  gage  too  snail 
covers  too  few  crystals  to  represent  the  behavior  of  the  material  on 
a  macroscopic  scale.  The  optimum  size  was  ascertained  to  be  a  J"  gage 
length,  and  the  gage  employed  was  a  Baldwin-Lima  Hamilton  paper-backed 
SR-4  type  FAF  transducer  with  a  resistance  of  120  ohms  and  a  gage 
factor  of  2.01.  These  gages  were  attached  on  opposite  ends  of  a 
diameter  and  coupled  in  series  to  eliminate  any  bending  components 
present.  The  gages  were  incorporated  in  an  AC-excited  bridge  circuit 
including  an  amplifier  and  the  output  was  recorded  on  a  plotter.  The 
records  were  calibrated  by  means  of  the  insertion  of  ten  known  shunt 
resistances  across  the  gage.  The  applied  compressive  force  was  deter¬ 
mined  by  means  of  a  calibrated  cell  with  a  0  -  10,000  lb  range  that 
concurrently  recorded  the  ordinate  of  the  stress-strain  curve. 

In  order  to  measure  the  presence  of  bending  which  shoulo  be  mini¬ 
mized  or  preferably  eliminated  from  these  static  tests,  the  two 
longitudinal  gages  were  decoupled  and  read  individually.  Initial  use 
of  various  types  of  end  pieces  including  a  bal  ~ocket,  rubber 

sheets  and  plaster  capping  ”ielded  deviations  3S  high  as  50  percent. 

The  best  results  were  obtained  with  the  aid  of  the  capping  material 
Hydrostone  which  consists  of  a  modified  Plaster  of  Paris  base.  This 
material  was  mixed  with  water  and  allowed  to  set  on  the  specimen  at  a 
pressure  of  20  psi  for  45  minutes.  This  reduced  the  f.eviation  between 
the  gage  readings  to  the  10-15  percent  range;  some  of  this  is  unavoidable 
due  to  local  inhomogeneities  and  the  nonrandom  orientation  of  crystals 
under  each  gage. 
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A  typical  compressive  quasi-static  stress-strain  curve  is 
presented  in  Figure  26,  clearly  exhibiting  the  nonlinear  response 
of  the  aaterial.  Figure  27  shows  the  first  ar.d  the  tenth  strain  cycle 
when  a  specimen  was  repeatedly  and  continuously  loaded  to  750  lb  and 
then  unloaded,  all  at  the  same  strain  rate.  Both  curves  indicate  the 
presence  of  hysteresis,  but  the  repeated  loading  apparently  diminishes 
this  effect  while  increasing  the  apparent  stiffness.  Further  tests 
along  this  line  under  varying  strain  rates  are  in  progress. 

The  development  of  tests  to  produce  intermediate  strain  rates  con¬ 
sisted  of  two  approaches,  the  first  of  which  attempted  to  transform  a 
long  rectangular  loading  pulse  produced  by  the  impact  of  a  long  cylin¬ 
der  on  a  Hopkinson  bar  into  a  triangular  pulse  which  exhibits  a 
predetermined  constant  strain  rate.  This  concept  employs  the  idea  of  a 
mechanical  converter  which  exploits  the  impedance  mismatch  of  concen¬ 
tric  bars  of  various  radii  and  material  properties  as  shown  in  Figure 
27  to  transform  such  a  rectangular  input  pulse  as  shown  in  Figure  28. 


The  transmission  and  reflection  coefficients,  T  and  R,  respectively, 

(12) 

for  this  converter  are  given  by 
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where  A  and  are  tne  cross-sectional  area  and  the  rod  wave  velocity 
of  the  various  elements,  respectively.  The  output  of  this  converter 
to  an  input  history  f(t)  is  given  by 


R(t)  = 


T  T 
12  23 
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n=0 


(R212  R232) 


ff  t 


~2 1 
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For  a  rectangular  input  waveform  of  unit  amplitude  and  duration  T* , 
choosing  £  to  be  such  that  cQ  =  I0i/T* ,  and  with  (R212R232^=  8 ’ 

T.„T  ,  =  0.035,  the  pulse  is  converted  into  an  approximately  triangular 

l<w  o 

shape  as  shown  in  Figure  29.  Other  input  histories  would  be  deformed 
in  an  analogous  fashion. 

Or.e  of  the  difficulties  associated  with  this  process  is  evident 
by  inspection  of  Figure  29,  namely  the  high  signal  loss  through  the 
converter.  It  may  also  be  noted  that  the  output  is  not  truly  ramp¬ 
shaped,  but  rather  multistepped  with  an  average  slope  which  decays  with 
t  ime. 


The  second  method  of  obtaining  intermediate  strain  rates  considered 
is  the  impact  of  a  long  cylinder  prefaced  by  a  cone  on  a  Hopkinson 
bar.  Numerous  materials  and  cone  geometries  of  the  composite  strik.r 
were  tested,  the  most  successful  results  being  obtained  with  a  30° 
steel  cone  mounted  on  the  end  of  a  J"  diameter,  8"  long  7')75 
aluminum  bar.  Strain  rates  of  10  sec  * ,  linear  ovi r  uic  entire  range 
of  strains  were  obtained  in  a  7075  aluminum  Hopkinson  bar,  as  shown 
in  Figure  30a.  In  a  standard  Hopkinson  bar  test,  the  strnn  rales  are 
of  tiie  order  of  1000  sec  *.  With  the  present  device,  the  rate  can  be 
reduced  from  this  value  only  to  a  level  of  about  50  see  1 ;  further 


ar 


reductions  require  the  use  of  a  different  technique.  The  use  of 
longer  input  bars  with  the  same  conical  striker  produced  an  undesirable 
distortion  of  pulse  from  its  triangular  shape  as  shown  in  Figure  30b. 
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III,  FUTURE  WORK 

Th  work  during  the  interval  up  to  February,  1973  will  be  devoted 
to  completing  the  major  objectives  of  the  program,  with  certain 
aspects  of  the  project  requiring  some  additional  time,  as  initially 
indicated  in  the  original  proposal.  Furthermore,  it  has  become  obvious 
that  the  techniques  developed  will  find  useful  applications  in  a  broader 
range  of  problems  in  the  field  of  rock  mechanics  than  the  scope  of 
the  work  covered  by  the  present  proposal.  Suggestions  for  such 
activities  as  a  sequel  to  the  current  investigation  will  be  present'd 
in  the  future. 

In  the  experimental  examination  of  wave  processes,  the  fo' lowing 
endeavors  are  planned:  (1)  The  development  an!  testing  of  new  and 
improved  grouting  materials  embedded  with  a  transducer  package  in  a 
Hopkinson  bar,  and  calibration  of  the  unit  to  determine  whether  the 
discontinuities  produced  by  drilling  have  been  satisfactorily  eliminated, 
(2)  The  drilling  of  deep  holes  in  Hopkinson  bars  of  Yule  marble  to 
simulate  the  installation  and  testing  of  the  transducer  packages  in 
the  rot k  blocks.  (3)  The  improvement  of  the  sensitivity  of  the  trans- 
'1m  i  r  puk.ige  by  the  use  of  ceramic  IZT4  crystals  that  exhibit  a 
hundredfold  larger  response  than  the  quartz  crystals  presently  employed. 
(1)  flic  employment  of  blocks  of  substitute  matt  rials  for  simulating  the 
drilling  and  embedment  process  of  th.-  transducers  in  the  actual  Yule 
marble  ha  1 1  -spaces ..  1 5 )  An  extensive  study  oi  the  Yule  marble  block 

geometry  to  locate  the  position  of  the  interior  t  r.insdi  ce  rs  and  their 
associated  drilling  holes  in  a  manner  mi  as  to  providt  tin  least  chanct 
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of  undeslred  mechanical  damage  to  the  block  and,  further,  to  provide 
maximum  utlity  for  the  gathering  of  ♦he  data  and  comparison  with 
theoretical  predictions.  (6)  The  investigation  of  the  use  of  trans¬ 
mission  devices  to  guard  against  damage  to  the  fabric  of  the  block  at 
the  impact  point.  (7)  Appropriate  checks  to  establish  the  repro¬ 
ducibility  of  the  data,  and  (8)  Improvement  of  existing  methods  of 
surface  measurement.  Future  investigations  may  even  involve  studies  of 
the  reflection  processes  from  bounding  surfaces  and  the  effect  of  other 
types  of  anisotropy. 

The  transform  analysis  will  be  concerned  with  the  further  verifi¬ 
cation  of  the  displacement  field  program,  the  development  of  more 
efficient  routines  for  extracting  roots  oi  the  characteristic  equation 
to  save  compute*"  time,  the  specification  of  a  convolution  integral 
procedure  so  as  to  permit  utilization  of  the  displacement  field  rela¬ 
tions  for  an  arbitrary  forcing  function,  and  the  development  of  a 
separate  program  for  the  evaluation  of  the  surface  response  which  at 
present  cannot  be  obtained  as  a  limiting  process  for  the  body  solution. 

In  order  to  permit  a  direct  comparison  with  experimental  results,  the 
response  of  the  block  must  be  predicted  in  terms  of  both  local  stress 
and  strain,  since  the  transducers  employed  m  the  experimental  phase 
of  the  program  measure  these  quantities.  This  requires  ,a  differentiation 
of  the  displacement  field  which  may  lead  to  additional  numerical  compli¬ 
cations  that  will  require  resolution.  It  will  be  necessary  to  closely 
correlate  the  computer  efforts  with  the  experimental  approach  so  as  to 
permit  the  calculation  of  results  at  positions  that  are  b >th  of  physical 
interest  and  experimentally  iecessibli,,  file  computation  can  be 
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considered  to  be  valid  only  at  a  given  station  for  time  intervals  prior 
to  the  arrival  of  significant  reflections  from  bounding  surfaces. 

The  finite  element  procedure  will  initially  concentrate  on  a 
better  representation  of  the  element  and  studies  concerned  with  the 
determination  of  the  optimal  spatial  mesh  density  and  time  step  size, 
using  an  isotropic  medium.  Upon  satisfactory  performance  of  the  pro¬ 
gram  in  the  testing  phases,  it  will  be  applied  to  the  subject  problem 
in  its  most  general  form.;  Reflections  from  the  boundary  will  also  be 
considered  in  additional  studies.  It  may  be  desirable  to  specialize 
this  program  for  the  case  of  plane  strain  in  a  special  anisotropic 

(7) 

medium  where  a  closed-form  solution  is  available  for  check  purposes 
The  numerical  results  obtained  here  be  compared  with  the 

corresponding  information  from  the  integral  transform  analysis. 

The  crystallographic  techniques  for  the  location  of  the  transversely 
isotropic  axis  of  the  Yule  marble  materials  will  be  further  pursued  and 
applied.  Tensile  and  compressive  tests  under  controlled  dynamic  condi¬ 
tions  will  be  executed  utilizing  the  concept  of  the  mechanical  converter 
previously  described,  and  standard  Hopkinson  bar  tests  on  the  material 
will  also  be  conducted  on  specimens  parallel  and  perpendicular  to  the 
axis  of  elastic  symmetry.:  A  major  effort  will  be  required  to  initiate 
the  program  of  fracture  analysis  of  the  material  as  a  function  of 

strain  rate.  For  this  purpose,  the  techniques  developed  on  bone, 

( 1  <>) 

utilizing  a  split  Hopkinson  bar,  will  be  employed  and  direct  observa¬ 

tion  of  the  fracture  process  by  means  of  high-^peed  ph ■ d og mphv  will  be 
sought.  The  tracing  of  the  crack  propagation  ni<  ms  iiansmitted 
light  effects  will  also  be  attempted. 


:9k 


49 


The  ultimate  goal  of  tins  portion  of  the  investigation  will  be 
the  complete  characterization  of  the  mechanical  response  of  this 
material  to  loads  ranging  in  rate  from  the  quasi-static  to  the  impact 
domain.  This  will  hopefully  be  accomplished  by  a  model  that  incorporates 
the  major  features  observed  during  the  current  experimental  program. 
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